#These three commands only need to be run once to install teh gpart packages
#chooseCRANmirror()
#install.packages("BiocManager")
#BiocManager::install("gpart")
library(gpart)

#Tell R that it can take parameters from the command line
args = commandArgs(trailingOnly=TRUE)
filename=args[1]
maxHap=args[2]

#Make some filenames from the name of the gene that was entered
gf=paste(filename, "ped",sep=".")
sf = paste(filename, "map",sep=".")
gif=paste(filename, "geneinfo",sep=".")
fn =paste(filename,"heatmap",sep=".")
outfile=paste(filename,"bigLD.txt",sep=".")

#Run BigLD to find haplotype blocks
res1 = BigLD(genofile = gf, SNPinfofile = sf, CLQcut = 0.5) 
write.table(res1, file = outfile, quote = FALSE, sep = "\t")

#GPART is good for larger regions but it did not work for me with single genes
#Gres0 = GPART(genofile = gf, SNPinfofile = sf, 
#              BigLDresult = res1, 
#              minsize = 4, maxsize = as.numeric(maxHap),
#              geneinfofile =gif)
              #geneDB = "ensembl", assembly = "GRCh37")

#You can generate a haplotype block image of the region
LDblockHeatmap(genofile = gf, SNPinfofile = sf,
               geneinfofile = gif, blockresult = res1,  
               filename = fn, minsize = 3, type = "tif",  res = 300)
